Ch/psestimate memory #91

cahofer · 2020-05-06T19:12:50Z

Monte Carlo Fisher matrix estimation for m-mode formalism for large telescope classes.

fix(psestimation): refactering the functions generate() and q_estimator() such they can be reused in psmc.MonteCarloXLarge.
feat(psmc.MonteCarloXLarge):
a class that estimates the Fisher matrix with Monte Carlo simulations when memory demand is high. Features:
- clarray for power spectrum bands is now an MPIarray distributed over bands. Fixes memory
  issues.
- the main function generate() divides all m-modes into chunks, so that each parallel process
  receives one m-mode at at time.
- the function recv_send_data() manages sending and receiving of m-modes from rank to
  rank. As such, each rank with a sub range of bands is subject to every m-mode.
- q_estimator() calculates the quadratic estimation under the assumption the q-array is a
  distributed MPIArray over bands

large telescopes

cahofer · 2020-07-09T19:14:23Z

Hey @jrs65! I added Tristan and Simon to the review. Simon is probably most familiar next to you with the Fisher estimation code, and Tristan is just a really good code reviewer. Let me know if there is someone else I could add to help with this. I would like to get this going and merged!

sjforeman

This looks good to me! (However, I'm not an expert code reviewer like @jrs65 and @tristpinsm , so it might be a good idea to have one of them skim through it before the final merge.)

My comments are all minor, and should be easily resolvable.

drift/core/psmc.py

drift/core/psestimation.py

drift/core/psmc.py

sjforeman · 2020-07-12T19:48:20Z

drift/core/psmc.py

+        return recv_buffer
+
+    def split_single(self, m, size):
+


Docstring would be useful here, to describe in words what this function does

Done. Need feedback if it is comprehensible.

Why use this routine instead of mpiutil.partition_list()?

I want every rank to work with one m-mode, with partition list it's not guaranteed that i-th partition will have be length = size, where size is the number of total MPI processes running.

Hmm, I think you could accomplish this with

mpiutil.partition_list( np.range(self.telescope.mmax + 1), mpiutil.rank, (self.telescope.mmax + 1) // size )

or something similar. However, if your code works, I'm not sure if it's worth the effort to change it...

jrs65

This is looking good, I'm glad to be finally getting this one merged.

I left a few comments, I'll have a look again when the current round of Simon's suggestions and my own are addressed. Does this actually run though? It seems like there's a few things which just shouldn't actually work.

drift/core/psmc.py

drift/core/psestimation.py

jrs65 · 2020-07-14T00:48:41Z

drift/core/psestimation.py

        else:
            y0 = x0
            y2 = x2

+        return x2, y2
+
+    def q_estimator(self, mi, vec1, vec2=None, noise=False):


Given that you have an entirely new implementation of q_estimator in your new class, why does the old one need to get changed?

I wanted to factor out the part of the code that projects from KL-basis to sky basis so that it can be reused as a method in the future.

and the new PS-code would get just like a big spaghetti code. Tried to factor out as much as I could.

Yep, makes sense to me

cahofer · 2020-12-05T04:05:59Z

This is looking good, I'm glad to be finally getting this one merged.

I left a few comments, I'll have a look again when the current round of Simon's suggestions and my own are addressed. Does this actually run though? It seems like there's a few things which just shouldn't actually work.

Yeah it's been running for half a year now...

cahofer · 2020-12-05T04:07:29Z

I did the changes you requested. @sjforeman @jrs65. I will run the code now and blacken before I commit and push the branch to git.

comments

lgtm-com · 2021-02-02T02:22:12Z

This pull request introduces 3 alerts and fixes 1 when merging 9f71ec6 into 94f0ece - view on LGTM.com

new alerts:

1 for Unused local variable
1 for Unused import
1 for Variable defined multiple times

fixed alerts:

1 for Unused local variable

cahofer · 2021-02-02T03:14:09Z

drift/core/psmc.py

+                if mi < (self.telescope.mmax + 1):
+                    self.q_estimator(mi, vec1)
+
+                    # TO DO: Calculate q_a for noise power (x0^H N x0 = |x0|^2)


I added her a TODO item - at the moment we have no capability of simulating a noise bias.
On a note:

the original psestimation code doesn't estimate a noise bias for a auto power spectrum (why?)

Is x0^2 (=normalized KL-modes) supposed to be the noise bias? According to the m-mode paper see equation (91) - I thought you need to sandwich the estimator between realizations of white noise to get b_a.

oh yeah, and the only that noise is set to True is in the Crosspower class and there it should be zero, so I don't see the point why it's even calculated?!

Comment to self and anyone else reading: We won't save noise in KL-basis to q-estimator. IMO it really doesn't need to be there.

I'm confused by these comments. Isn't the noise bias estimated by the combination of fisher_bias and gen_sample?

I guess the problem is that the noiseonly argument isn't set to True when zero_mean == False which it should be?

cahofer · 2021-02-02T03:16:57Z

drift/core/psmc.py

+                        axis=0,
+                    ).astype(np.float64)
+
+    def noise_bias(self, mi, vec1, vec2=None):


I added this method for the future, in case we ever want to have a routine to calculate what Richard called the noise bias.

Scrap that. I don't think we need it.

cahofer · 2021-02-02T03:18:52Z

drift/core/psmc.py

+            request.Wait()
+
+            # Fill recv buffer with messages
+            recv_buffer[slc] = bi


oh yeah, I have been getting a warning that this should be recv_buffer[tuple(slc)] instead of recv_buffer[slc]. Any opinions?

cahofer

Okay, this is the most up to date PR for the multi parallelized Fisher code. It is ready for a second round. I left some comments concerning calculating a noise bias.
The new code as it is works at the moment only for

auto correlation Fisher errors
when no noise bias is calculated (but neither did the old psestimation code.)

cahofer · 2021-02-02T03:24:42Z

Oh - and for some reason the black check is failing right now even though I reformatted the files with black - so I am not sure whats going on there.

lgtm-com · 2021-02-02T03:24:53Z

This pull request fixes 2 alerts when merging c72303e into 94f0ece - view on LGTM.com

fixed alerts:

2 for Unused local variable

sjforeman

I have a few minor comments, but I will punt some of the more technical points to @jrs65

drift/core/psestimation.py

drift/core/beamtransfer.py

sjforeman · 2021-02-11T22:16:37Z

drift/core/psestimation.py

        else:
            y0 = x0
            y2 = x2

+        return x2, y2
+
+    def q_estimator(self, mi, vec1, vec2=None, noise=False):


Yep, makes sense to me

drift/core/psestimation.py

sjforeman · 2021-02-11T22:23:52Z

drift/core/psmc.py

@@ -81,7 +85,7 @@ def _work_fisher_bias_m(self, mi):
            x = self.gen_sample(mi, n)
            qa[:, s:e] = self.q_estimator(mi, x)

-        ft = np.cov(qa)
+        # ft = np.cov(qa)


Looks like this could be removed completely

I agree, I will get rid of it... I just kept it because it was Richard's comments.

drift/core/psmc.py

sjforeman · 2021-02-11T22:42:52Z

drift/core/psmc.py

+        return recv_buffer
+
+    def split_single(self, m, size):
+


Why use this routine instead of mpiutil.partition_list()?

jrs65 · 2021-05-18T02:27:49Z

@cahofer it looks like there's a few more comments being worked on in here, is that right?

sjforeman · 2021-06-02T01:07:50Z

drift/core/psestimation.py

-    def q_estimator(self, mi, vec1, vec2=None, noise=False):
-        """Estimate the q-parameters from given data (see paper).
-
+    def project_vector_kl_to_sky(self, mi, vec1):


Same comment as before: should this go in drift.core.kltransform instead of drift.core.psestimation?

…y/driftscan into ch/psestimate-memory

Carolin Hofer added 7 commits April 1, 2020 15:36

fix(psestimation): refactor generate and q_estimator

9ab5353

feat(PSMonteCarloXLarge): Fisher matrix via Monte-Carlo simulations for

5b88074

large telescopes

bug fixes

05ae3d4

blacken and some more bug fices

091a1d6

blacken

dad3e71

fix(psmc): cleaned up MonteCarloXlarge

3580934

blacken

6e6ad1e

cahofer requested a review from jrs65 May 6, 2020 19:12

cahofer mentioned this pull request May 6, 2020

Ch/psestimate memory radiocosmology/draco#96

Open

cahofer requested review from sjforeman and tristpinsm July 9, 2020 19:12

sjforeman requested changes Jul 12, 2020

View reviewed changes

jrs65 requested changes Jul 14, 2020

View reviewed changes

drift/core/psmc.py Outdated Show resolved Hide resolved

drift/core/psmc.py Outdated Show resolved Hide resolved

drift/core/psmc.py Outdated Show resolved Hide resolved

drift/core/psestimation.py Outdated Show resolved Hide resolved

jrs65 requested changes Jul 14, 2020

View reviewed changes

Carolin Hofer and others added 5 commits January 14, 2021 20:47

fix(psmc, psestimation, manager, kltransform): added docstrings and

eac383f

comments

docs(psestimation, psmc): update documentation in PSMonteCarloLarge

c466098

fix(psestimation,psmc): docstrings and variable names

c18e4fc

Merge branch 'master' into ch/psestimate-memory

2f49c17

fix(psmc): check if y2 is set to None in method q_estimator

9f71ec6

cahofer force-pushed the ch/psestimate-memory branch from 40435ab to 9f71ec6 Compare February 2, 2021 01:58

fix(psmc): black and fixing lgtm alerts

c72303e

cahofer force-pushed the ch/psestimate-memory branch from a9a76a2 to c72303e Compare February 2, 2021 03:01

cahofer commented Feb 2, 2021

View reviewed changes

sjforeman requested changes Feb 11, 2021

View reviewed changes

sjforeman reviewed Jun 2, 2021

View reviewed changes

Carolin Hofer and others added 3 commits June 19, 2021 14:48

fix : code cleanup, option for noise bias estimate

8940660

Merge branch 'ch/psestimate-memory' of ssh://github.com/radiocosmolog…

552493d

…y/driftscan into ch/psestimate-memory

Merge branch 'master' into ch/psestimate-memory

bc64134

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Ch/psestimate memory #91

Ch/psestimate memory #91

cahofer commented May 6, 2020

cahofer commented Jul 9, 2020

sjforeman left a comment

sjforeman Jul 12, 2020

cahofer Dec 5, 2020

sjforeman Feb 11, 2021

cahofer May 24, 2021

sjforeman Jun 2, 2021

jrs65 left a comment

jrs65 Jul 14, 2020

cahofer Dec 5, 2020

cahofer Feb 2, 2021

sjforeman Feb 11, 2021

cahofer commented Dec 5, 2020

cahofer commented Dec 5, 2020

lgtm-com bot commented Feb 2, 2021

cahofer Feb 2, 2021

cahofer Feb 2, 2021

cahofer May 24, 2021

jrs65 Jun 2, 2021

cahofer Feb 2, 2021

cahofer May 24, 2021

cahofer Feb 2, 2021

cahofer left a comment

cahofer commented Feb 2, 2021

lgtm-com bot commented Feb 2, 2021

sjforeman left a comment

sjforeman Feb 11, 2021

sjforeman Feb 11, 2021

cahofer May 24, 2021

sjforeman Feb 11, 2021

jrs65 commented May 18, 2021

sjforeman Jun 2, 2021

Ch/psestimate memory #91

Are you sure you want to change the base?

Ch/psestimate memory #91

Conversation

cahofer commented May 6, 2020

cahofer commented Jul 9, 2020

sjforeman left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jrs65 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

cahofer commented Dec 5, 2020

cahofer commented Dec 5, 2020

lgtm-com bot commented Feb 2, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

cahofer left a comment

Choose a reason for hiding this comment

cahofer commented Feb 2, 2021

lgtm-com bot commented Feb 2, 2021

sjforeman left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jrs65 commented May 18, 2021

Choose a reason for hiding this comment